Setup and Read Data

source("ams_initialize_script.R")
source("SCIM_calculation.R")  
source("ivsc_2cmt_RR_V1.R")
dirs$rscript_name = "Task16b_Parallel_Coordinates_Soluble_2019-11-12_AFIRvsSCIM_T0assumpt.R"
dirs$filename_prefix= str_extract(dirs$rscript_name,"^Task\\d\\d\\w?_")

data_in = read.csv("results/Task15_2019-11-12_10e3.csv",stringsAsFactors = FALSE)
data_in$id = 1:1e4

Compute various quantities for comparing AFIR and SCIM, theory and simulation

Using adhoc theory calculation for SCIM

#put data into categories ----
data    = data_in %>%
  mutate(SCIM_thy= SCIM_adhoc_thy) %>%
  mutate(Cavgss  = dose_nmol/(CL*tau),
         Kss_TL  = Kd_TL + keTL/kon_TL,
         Kss_DT  = Kd_DT + keDT/kon_DT,
         TL0     = T0*L0/Kss_TL,
         koff_TL = Kd_TL*kon_TL,
         koff_DT = Kd_DT*kon_DT,
         ksynT   = T0*keT + keTL*TL0,
         ksynL   = L0*(kon_TL*T0 + keL) - koff_TL*TL0,
         Lss     = ksynL/keL) %>%
  mutate(AFIR_SCIM_sqerr = (AFIR_thy - SCIM_sim)^2) %>%
  mutate(AFIRthy_category = case_when(AFIR_thy <  0.05 ~ "AFIRthy < 5%",
                                      AFIR_thy >  0.30 ~ "AFIRthy > 30%",
                                      AFIR_thy >= 0.05 &  AFIR_thy <= 0.30 ~ "5% <= AFIRthy <= 30%"),
         AFIRsim_category = case_when(AFIR_sim <  0.05 ~ "AFIRsim < 5%",
                                      AFIR_sim >  0.30 ~ "AFIRsim > 30%",
                                      AFIR_sim >= 0.05 &  AFIR_sim  <= 0.30 ~ "5% <= AFIRsim <= 30%"),
         SCIMthy_category = case_when(SCIM_thy <  0.05 ~ "SCIMthy < 5%",
                                      SCIM_thy >  0.30 ~ "SCIMthy > 30%",
                                      SCIM_thy >= 0.05 &  SCIM_thy  <= 0.30 ~ "5% <= SCIMthy <= 30%"),
         SCIMsim_category = case_when(SCIM_sim <  0.05 ~ "SCIMsim < 5%",
                                      SCIM_sim >  0.30 ~ "SCIMsim > 30%",
                                      SCIM_sim >= 0.05 &  SCIM_sim  <= 0.30 ~ "5% <= SCIMsim <= 30%"),
         AFIRthy_AFIRsim_category = paste0(AFIRthy_category, ", ", AFIRsim_category),
         AFIRthy_SCIMsim_category = paste0(AFIRthy_category, ", ", SCIMsim_category),
         AFIRsim_SCIMsim_category = paste0(AFIRsim_category, ", ", SCIMsim_category),
         SCIMthy_SCIMsim_category = paste0(SCIMthy_category, ", ", SCIMsim_category),
         error_category = case_when(AFIR_SCIM_sqerr < 0.1 ~ "low_error",
                                    TRUE ~ "high_error"))
data = data %>%
  arrange(AFIR_thy) %>%
  mutate(AFIRthy_category = factor(AFIRthy_category, levels = unique(AFIRthy_category))) %>%
  arrange(AFIR_sim) %>%
  mutate(AFIRsim_category = factor(AFIRsim_category, levels = unique(AFIRsim_category))) %>%
  arrange(SCIM_sim) %>%
  mutate(SCIMsim_category = factor(SCIMsim_category, levels = unique(SCIMsim_category)))

#check the assumptions of the data ----
data = data %>%
  mutate(Ttotss                 = T0*Tfold,
         koff_DT                = Kd_DT*kon_DT,
         assumption_plateau     = AFIR_thy < 0.30,
         assumption_drug_gg_T0  = Cavgss > 10*Ttotss,
         assumption_drug_gg_KssDT = Cavgss > 10*Kss_DT,
         assumption_koffDT_gt_keT = koff_DT > keT,
         assumption_koffTL_fast   = koff_TL > 1/30,
         assumption_Cavgss_gg_LssKssDT_KssTL = Cavgss > 10*Kss_DT*Lss/Kss_TL,
         assumption_T0simple    = T0/(ksynT/keT) > 0.5 & T0/(ksynT/keT) < 2, #the simple formula works for T0
         assumption_Lconst      = Lss/L0 > 0.5 & L0/Lss < 2, #then SCIM = AFIR
         assumption_Tss_gt_Lss  = Tss_sim > Lss_sim,
         assumption_all         = assumption_plateau & 
                                  assumption_drug_gg_T0 &
                                  assumption_drug_gg_KssDT &
                                  assumption_koffDT_gt_keT & 
                                  assumption_koffTL_fast &           
                                  assumption_Cavgss_gg_LssKssDT_KssTL &
                                  assumption_T0simple &
                                  assumption_Tss_gt_Lss &
                                  assumption_Lconst,
         assumption_category    = "",
         assumption_category    = ifelse(assumption_plateau,"","!plateau, "),
         assumption_category    = ifelse(assumption_drug_gg_T0,"","drug !>> T0"))

Put data into error categories and summarize

threshold = 0.1
data_errss = data_in %>%
  filter(abs(TLss_frac_change)>=threshold) 
print(paste0(nrow(data_errss)," of ", nrow(data_in), " : Number of rows with TLss_frac_change > 0.1"))
## [1] "208 of 10000 : Number of rows with TLss_frac_change > 0.1"
data_err0 = data_in %>%
  filter(abs(TL0_05tau_frac_change)>=threshold)
print(paste0(nrow(data_err0)," of ", nrow(data_in), " : Number of rows with TL0_05tau_frac_change > 0.1"))
## [1] "0 of 10000 : Number of rows with TL0_05tau_frac_change > 0.1"
# error historgram ----
data_quick_summ = data %>%
  select(id,AFIR_thy, SCIM_sim, AFIR_SCIM_sqerr, TLss_frac_change, TL0_05tau_frac_change) %>%
  gather(key,value,-c(id)) %>%
  mutate(category = case_when((value < threshold) ~ "keep_low",
                              ((value >= threshold) & (key %in% c("AFIR_SCIM_sqerr","SCIM_sim"))) ~ "keep_high",
                              ((value >= threshold) & (key %in% c("AFIR_thy"))) ~ "keep_high_AFIR",
                              TRUE ~ "remove_high_error"))

g = ggplot(data_quick_summ, aes(value, fill = category))
g = g + geom_histogram()
g = g + facet_wrap(~key, scales = "free")
g = g + scale_fill_manual(values = c(keep_low = "grey80", 
                                     keep_high = "grey50",
                                     remove_high_error = "red", 
                                     keep_high_AFIR = "blue"))
g = g + xgx_scale_x_log10()
g = g + ggtitle("")
print(g)

#keep only the simulations with no issues 
data_keep = data %>%
  filter(TLss_frac_change < threshold, 
         TL0_05tau_frac_change < threshold)

#put simulations into different categories
data_summary = data_keep %>%
  group_by(AFIRthy_SCIMsim_category) %>%
  count() %>%
  arrange(desc(n))
kable(data_summary)
AFIRthy_SCIMsim_category n
AFIRthy < 5%, SCIMsim < 5% 4270
AFIRthy < 5%, SCIMsim > 30% 1667
5% <= AFIRthy <= 30%, SCIMsim > 30% 1275
5% <= AFIRthy <= 30%, 5% <= SCIMsim <= 30% 997
AFIRthy > 30%, SCIMsim > 30% 916
AFIRthy < 5%, 5% <= SCIMsim <= 30% 759

AFIRsim vs SCIMsim : 3x3 plot colors —-

param2uniform = function(x) {(log(x) - log(min(x)))/(log(max(x))-log(min(x)))}
data_plot = data_keep %>%
  mutate_at(vars(AFIR:kon_TL,dose_mpk), funs(tf=param2uniform(.))) %>%
  select(id,contains("AFIR"),contains("SCIM"), T0_tf:kon_TL_tf, dose_mpk_tf, contains("assumption")) %>%
  gather(param,param_value,-c(id, contains("AFIR"), contains("SCIM"), contains("assumption"))) %>%
  mutate(param = str_replace(param,"_tf",""))

#sort by average param value in one category to help with visualization ----
data_summ = data_plot %>%
  filter(AFIRthy_SCIMsim_category == "AFIRthy < 5%, SCIMsim > 30%") %>%
  group_by(param,AFIRthy_SCIMsim_category) %>%
  summarise(x = mean(param_value)) %>%
  arrange(x) %>%
  ungroup()
kable(data_summ)
param AFIRthy_SCIMsim_category x
Kd_TL AFIRthy < 5%, SCIMsim > 30% 0.3746478
keT AFIRthy < 5%, SCIMsim > 30% 0.4261922
Kd_DT AFIRthy < 5%, SCIMsim > 30% 0.4417119
Lfold AFIRthy < 5%, SCIMsim > 30% 0.4664135
dose_mpk AFIRthy < 5%, SCIMsim > 30% 0.4688410
kon_DT AFIRthy < 5%, SCIMsim > 30% 0.5014002
keL AFIRthy < 5%, SCIMsim > 30% 0.5076666
keDT AFIRthy < 5%, SCIMsim > 30% 0.5169649
kon_TL AFIRthy < 5%, SCIMsim > 30% 0.5849314
L0 AFIRthy < 5%, SCIMsim > 30% 0.6468130
T0 AFIRthy < 5%, SCIMsim > 30% 0.6970433
data_plot = data_plot %>%
  mutate(param = factor(param, 
                        levels = data_summ$param))

g = ggplot(data_plot, aes(x=param,y=param_value, group = id, color = assumption_all, alpha = assumption_all))
g = g + geom_line()
g = g + facet_grid(SCIMsim_category~AFIRthy_category,switch = "y")
g = g + theme(axis.text.x = element_text(angle = 45, hjust = 1))
g = g + labs(x = "Parameter", y = "Parameter Value")
g = g + guides(colour = guide_legend(override.aes = list(alpha = 1)))
g = g + scale_color_manual(values = c(`TRUE` = "blue", `FALSE` = "red"))
g = g + scale_alpha_manual(values = c(`TRUE` = .1, `FALSE` = 0.01))
g = xgx_save(7,7,dirs,"Parallel_Coord_Soluble_3x3_AFIRthy_AFIRsim","")
g1 = g
print(g)

Focus on when assumptions are true and SCIM > 30% while AFIRthy < 5%

#explore data data where all assumptions are true and still
#AFIRsim > 30% and AFIRthy < 5% ---- on look, there is lots of L0!!!
#focus on this plot
data_new = data_plot %>%
  filter(SCIMsim_category == "SCIMsim > 30%",
         AFIRthy_category == "AFIRthy < 5%",
         assumption_all == TRUE)
g = g1
g = g %+% data_new
g = g + geom_line(alpha = 0.05)
g = xgx_save(5,5,dirs,"Parallel_Coord_Soluble_AFIRthy_lt_5_SCIMsim_ge_30","")
g2= g
#print(g)

Identify patients where all assumptions true and theory vs sim disagree.

source("SCIM_calculation.R")  
source("ivsc_2cmt_RR_V1.R")
model = ivsc_2cmt_RR_KdT0L0()
library(RxODE)

id = unique(sort(data_new$id))
print("these IDs, even with all the restrictions, AFIR and SCIM still don't match")
## [1] "these IDs, even with all the restrictions, AFIR and SCIM still don't match"
print(id)
##  [1]    6  267 1289 2519 3186 4526 4621 5355 5484 5838 6009 6268 6322 7077
## [15] 7593 8379 8462 9268 9302 9624
for (id_plot in id[1:5]) {

  g = g2
  g = g + geom_line(data = filter(data_new,id==id_plot),
                    size = 2,
                    color = "black")
  g = g + ggtitle(paste("id =",id_plot))
  print(g)
  filepref = paste0("Parallel_Coord_Soluble_AFIRthy_lt_5_SCIMsim_ge_30_",id_plot)
  g = xgx_save(5,5,dirs,filepref,"")
  
  #simulate a patient where theory and simulation disagree 
  param = data %>%
    filter(id==id_plot)
  
  assumptions = param %>%
    select(contains("assumption")) %>%
    t() 
  kable(assumptions)
  
  tmax = 365 #days
  tau  = param$tau   #days
  dose_nmol = param$dose_nmol
  compartment = 2
  infusion = TRUE
  
  nam   = names(param)
  param_as_double = param %>%
    as.numeric() %>%
    setNames(nam)
  param_as_double = param_as_double[model$pin]
  
  param_print = param_as_double %>%
    t() %>%
    as.data.frame() %>%
    mutate(CL = signif(keD/V1,2),
           id = id_plot) %>%
    select(id, CL,T0,L0,Kd_DT,Kd_TL,kon_DT,kon_TL,keT,keL,keDT,keTL)
  
  ev = eventTable(amount.units="nmol", time.units="days")
  sample.points = c(seq(0, tmax, 0.1), 10^(-3:0)) # sample time, increment by 0.1
  sample.points = sort(sample.points)
  sample.points = unique(sample.points)
  ev$add.sampling(sample.points)
  if (infusion == FALSE) {
    ev$add.dosing(dose=dose_nmol, start.time = tau, nbr.doses=floor(tmax/tau), dosing.interval=tau, dosing.to=compartment)
  } else {
    ev$add.dosing(dose=dose_nmol, start.time = tau, nbr.doses=floor(tmax/tau)+1, dosing.interval=tau, dosing.to=compartment, dur = tau)
  }  
  
  sim = lumped.parameters.simulation(model, param_as_double, dose_nmol, tmax, tau, compartment, infusion)
  thy = lumped.parameters.theory    (       param_as_double, dose_nmol,       tau,              infusion)
  
  sim_rename = sim
  nam = names(sim_rename) %>%
    str_replace_all("_sim$","")
  names(sim_rename) = nam
  sim_rename$type = "sim"
  
  thy_rename = thy
  nam = names(thy_rename) %>%
    str_replace_all("_thy$","")
  names(thy_rename) = nam
  thy_rename$type = "thy"
  
  compare = bind_rows(sim_rename,thy_rename) %>%
    select(type,Dss,T0,L0,TL0,Ttotss,Lss,TLss,AFIR,SCIM)

  init = model$init(param_as_double)
  out  = model$rxode$solve(model$repar(param_as_double), ev, init)
  out  = model$rxout(out)
  
  out_plot = out %>%
    select(time,D,T,DT,L,TL) %>%
    gather(cmt,value,-time)
  out_last = out_plot[(out$time==max(out$time)),]
  
  g = ggplot(out_plot,aes(x=time,y=value, color = cmt, group= cmt))
  g = g + geom_line()
  g = g + geom_label(data = out_last, aes(label = cmt), show.legend = FALSE, hjust=1)
  g = g + geom_vline(xintercept = tau, linetype = "dotted")
  g = g + xgx_scale_x_time_units(units_dataset = "days", units_plot = "weeks")
  g = g + xgx_scale_y_log10()
  g = g + labs(y = "Concentration (nm)", color = "")
  g = g + ggtitle(paste0(  "id = ",param$id,
                           "\nAFIR_thy  = ",signif(thy$AFIR_thy,2),
                           "\nAFIR_sim  = ",signif(sim$AFIR_sim,2),
                           "\nSCIM_thy = ",signif(thy$SCIM_adhoc_thy,2),
                           "\nSCIM_sim = ",signif(sim$SCIM_sim,2)))
  filepref = paste0("Parallel_Coord_Soluble_AFIRthy_lt_5_SCIMsim_ge_30_",id_plot)
  g = xgx_save(5,5,dirs,filepref,"")
  print(g)
  
  #unfortunately, kable does not work properly inside for loop
  print(t(param_print))
  print(t(compare))
}

##            [,1]
## id     6.00e+00
## CL     2.70e-01
## T0     2.58e+01
## L0     1.57e-02
## Kd_DT  1.61e+01
## Kd_TL  1.04e-01
## kon_DT 3.83e+02
## kon_TL 2.82e+01
## keT    1.32e+03
## keL    4.40e+00
## keDT   7.06e-02
## keTL   1.27e+00
##        [,1]         [,2]        
## type   "sim"        "thy"       
## Dss    "6645781"    "6686508"   
## T0     "25.8"       "25.8"      
## L0     "0.0157"     "0.0157"    
## TL0    "2.717877"   "2.717877"  
## Ttotss "461494.0"   "482428.5"  
## Lss    "0.2528003"  "0.8001780" 
## TLss   "  1.896427" "-99.000000"
## AFIR   "0.04333395" "0.04308006"
## SCIM   "0.6977604"  "2.2946954"

##            [,1]
## id     2.67e+02
## CL     2.70e-01
## T0     1.61e+01
## L0     5.44e-03
## Kd_DT  8.58e+00
## Kd_TL  6.55e-01
## kon_DT 1.53e+03
## kon_TL 5.60e+00
## keT    1.17e+03
## keL    5.92e+00
## keDT   8.16e-02
## keTL   4.76e+00
##        [,1]           [,2]          
## type   "sim"          "thy"         
## Dss    "2457735"      "2480159"     
## T0     "16.1"         "16.1"        
## L0     "0.00544"      "0.00544"     
## TL0    "0.05819535"   "0.05819535"  
## Ttotss "219844.3"     "230849.0"    
## Lss    "0.03704329"   "0.05223221"  
## TLss   "  0.01889042" "-99.00000000"
## AFIR   "0.04766974"   "0.04725864"  
## SCIM   "0.3246037"    "0.4762670"

##             [,1]
## id     1.289e+03
## CL     2.700e-01
## T0     1.410e+02
## L0     9.100e-04
## Kd_DT  6.030e+01
## Kd_TL  1.040e-02
## kon_DT 2.160e+01
## kon_TL 3.500e+03
## keT    1.530e+00
## keL    3.130e+00
## keDT   4.520e-02
## keTL   3.770e-01
##        [,1]          [,2]         
## type   "sim"         "thy"        
## Dss    "319176.0"    "319444.4"   
## T0     "141"         "141"        
## L0     "0.00091"     "0.00091"    
## TL0    "12.21103"    "12.21103"   
## Ttotss "4763.398"    "4874.636"   
## Lss    "0.1303555"   "1.4716952"  
## TLss   " 11.13632"   "-99.00000"  
## AFIR   "0.006366509" "0.006349270"
## SCIM   " 0.9119888"  "10.2807213"

##             [,1]
## id     2.519e+03
## CL     2.700e-01
## T0     1.990e+01
## L0     9.150e-04
## Kd_DT  1.420e+01
## Kd_TL  2.160e-03
## kon_DT 5.690e+02
## kon_TL 6.360e+02
## keT    2.350e+02
## keL    4.170e+01
## keDT   1.020e-01
## keTL   4.480e+01
##        [,1]          [,2]         
## type   "sim"         "thy"        
## Dss    "5867197"     "5873016"    
## T0     "19.9"        "19.9"       
## L0     "0.000915"    "0.000915"   
## TL0    "0.2508049"   "0.2508049"  
## Ttotss "45635.42"    "45958.20"   
## Lss    "0.1026274"   "0.2703649"  
## TLss   "  0.1561306" "-99.0000000"
## AFIR   "0.005550212" "0.005539720"
## SCIM   "0.6225181"   "1.6499348"

##             [,1]
## id     3.186e+03
## CL     2.700e-01
## T0     5.680e+01
## L0     1.550e-03
## Kd_DT  7.720e+01
## Kd_TL  1.050e-01
## kon_DT 5.080e+00
## kon_TL 1.760e+01
## keT    2.540e+01
## keL    1.980e+00
## keDT   6.660e-02
## keTL   2.660e-01
##        [,1]          [,2]         
## type   "sim"         "thy"        
## Dss    "821671.6"    "823412.7"   
## T0     "56.8"        "56.8"       
## L0     "0.00155"     "0.00155"    
## TL0    "0.7329726"   "0.7329726"  
## Ttotss "20916.30"    "21665.39"   
## Lss    "0.03127469"  "0.10002005" 
## TLss   "  0.5117136" "-99.0000000"
## AFIR   "0.03460014"  "0.03452809" 
## SCIM   "0.6981347"   "2.3076631"